Nonlinear coherent transport of waves in disordered media 
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We present a diagrammatic theory for coherent backscattering from disordered dilute media in 
the nonlinear regime. The approach is non-perturbative in the strength of the nonlinearity. We 
show that the coherent backscattering enhancement factor is strongly affected by the nonlinearity, 
and corroborate these results by numerical simulations. Our theory can be applied to several 
physical scenarios like scattering of light in nonlinear Kerr media, or propagation of matter waves 
in disordered potentials. 
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It is already known that the interplay between disor- 
der and - even very weak - nonlinearity can lead to dra- 
matic changes to the system's properties: for example, 
instabilities occur P, 0, S], or localization may be de- 
stroyed Q. In the experiments studying the localization 
properties of matter waves in speckle potentials [51] i the 
nonlinear regime, arising from the atomic interactions, 
is almost unavoidable. Furthermore, nonlinear behavior 
is easily observed in coherent backscattering experiments 
using cold atoms as scatterers 6ij. As a third example, 
also the random laser exhibits nonlinearities which po- 
tentially influence the structure of localized laser modes 
0. In ah these cases, even if the systems are governed by 
simple nonlinear wave equations, a precise description of 
the impact of this nonlinearity on the interference effects 
altering the properties of diffuse wave propagation is still 
lacking. Since exact numerical calculations for realistic 
situations are at the border of or beyond actual com- 
puter capacities, one needs an efhcient theory providing 
directly disorder averaged quantities. For this purpose, 
the present letter shows that the standard diagrammatic 
approach [1, [^, can be extended to the nonlinear 
regime. Using ladder and crossed-like diagrams, we will 
derive a nonlinear radiative transfer equation for the av- 
eraged wave intensity, and then calculate the interference 
corrections on top of the nonlinear solution. 

The general frame where our approach can be applied 
is as follows: we assume a nonlinear wave equation with 
unique and stationary monochromatic solution. In par- 
ticular, we assume that all even orders x^"''? n — 2, 4, ... , 
of the nonlinear susceptibility vanish, such that the gen- 
eration of higher harmonics can be neglected. Further- 
more, the refractive index modifications are small enough 
such that we can neglect effects like self-focusing, pattern 
formation and solitons [ll| on the length scale set by 
the disorder (a mean free path). Instead, the nonlinear 
effects relevant for our disordered case can be summa- 
rized as follows: firstly, the wave intensity /(r) becomes a 
fluctuating quantity, which is especially important in the 
nonlinear regime. Secondly, the usual picture of weak lo- 
calization resulting from interference only between pairs 
of amplitudes propagating along reversed paths breaks 



down in the nonlinear regime. As a consequence of the 
nonlinear mixing between different partial waves, weak 
localization must rather be interpreted as a multi-wave 
interference phenomenon [12,. d3.]. In particular, we will 
show in the following that the height of the coherent 
backscattering peak is strongly affected by nonlinearities, 
even if they do respect the reci pro city symmetry. In con- 
trast to our previous work (l^. Il3l| , the present approach 
is non-perturbative in the strength of the nonlinearity. 

At first, we consider an assembly of N point-like scat- 
terers located at randomly chosen positions r^, « = 
1, . . . , inside a sample volume V illuminated by a plane 
wave ki, with jkij = k. We assume the field radiated 
by each scatterer to be a nonlinear function f{Ei) of the 
local field Ei. Since all even orders of the nonlinear sus- 
ceptibilities vanish, we can write f{E) = g{I)E, where 
/ = EE* is the local intensity, and g{I) is proportional 
to the polarizability of the scatterers. This results in a 
set of nonlinear equations for the field at each scatterer: 
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As explained above, we aim at providing a theory pro- 
viding the relevant quantities (local intensities, coherent 
backscattering cone...) averaged over the random posi- 
tions of the scatterers. In a first step, we will derive 
an equation for the mean intensity (/(r)). In the dilute 
regime, where the typical distances jr^ — rj\ are much 
larger than the wavelength, we may neglect - in first ap- 
proximation - correlations between the fields emitted by 
different scatterers. Under this condition, the scattered 
field Ed(v) is a superposition of spherical waves with ran- 
dom relative phases, depicting thus a speckle pattern. 
This leads to the well known Gaussian statistics for the 
complex field Ed{r) [13], which is thus completely deter- 
mined by a single parameter, the mean diffuse intensity 
/d(r) = {\Ed{r)\'^). In addition to the scattered field, 
there is also a non-fluctuating coherent component origi- 
nating directly from the incident field. In total, we have 
E{r) — Ec{r) + Ed{r), and the average intensity sphts 
into a coherent and diffuse part: (/(r)) = /c(r) + /rf(r). 
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with Ic = I -Be P • The mean density of radiation intensity 
emitted from point r is then given by: 



i^(r)=AA(/r)=AA(|.g(/(r))p/(r)) 



(2) 



where M = N /V denotes the density of scatterers, and 
the average (...) is taken over the Gaussian statistics of 
the scattered field. 

Between two scattering events, the wave propagates in 
an effective medium made by the scatterers. The corre- 
sponding refractive index n and mean free path ^ are not 
the same for the coherent and the diffuse fields, because 
of their different statistical properties combined with the 
nonlinear behavior of the scatterers. In the dilute regime, 
the diffuse amplitude can be considered as a weak probe, 
such that the complex refraction index reads: 



n = l 
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whereas, for the coherent mode, the derivative d/ dE is 
replaced by l/£'c, i.e. = 1 + (f) / {2k'^ E^) , and 
l/£c = 2fclm{nc}. Since the results of the averages de- 
pend on /c(r) and /d(r), the nonlinear refractive indices 
also attain a spatial dependence n(r) and ndr). They 
describe average propagation of one strong and many un- 
correlated weak fields. (If there is more than one strong 
field, additional phenomena like four-wave mixing occur 

Recollecting all preceding ingredients, the transport 
equations for the average intensity read as follows: 
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Here, z denotes the distance from the surface of V 
to r, in the direction of the incident beam. Further- 
more, propagation from r' to r implies a spatial aver- 
age of l/£{r), which we note as follows: \r — r'\/i := 
r — r'l /p ds/£{r--sr + sr'), and similarly for the coherent 

mode {z/£c)- Since K, £ and £c depend on /c(r) and /(i(r), 
the above Eqs. (I4|5p form two coupled integral equations, 
whose solution we find numerically. Finally, the intensity 
scattered into backwards direction (expressed as dimen- 
sionless quantity, the so-called 'bistatic coefficient' re- 
sults as: 
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where A denotes the transverse (with respect to the in- 
cident beam) area of the scattering volume V. 

The validity of the preceding approach has been tested 
using the following nonlinear function: 
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Figure 1: (Color online) Comparison between exact numer- 
ical calculations and the theoretical approach (see text for 
details). The solid line depict the exact numerical results, 
whereas the dashed line corresponds to Fl including geomet- 
rical effects. The dotted line corresponds to the sum Fl -l-Fc 
exactly in the backward direction. The additional curve (long 
dashed line) plotted for a = 0.2 depicts the results obtained 
when the fluctuating character of the diffuse field is not taken 
into account. 



which depicts the (elastic) nonlinear behavior of a two- 
level atom exposed to an intense laser beam. The non- 
linear scatterers described by Eq. ([7]) are randomly dis- 
tributed inside a sphere, with homogeneous density. We 
must emphasize that, for this particular model of non- 
linearity, the stationary solution is always found to be 
unique and stable, as a consequence of the saturation 
g{I) —^ for large a. From the numerical solution of 
Eq. Il]), we calculate the radiated field and intensity out- 
side the cloud in different directions 9. This procedure is 
then repeated with many different configurations giving 
us the disorder averaged field and intensity. The results 
presented in this letter are obtained with 3000 configu- 
rations of 1500 scatterers with density such that k£ — 67 
and optical thickness 6 = 2 (in the linear limit a = 0). 

The results for the average intensity as a function of 
the backscattering angle 9 are depicted in Fig. [T] for differ- 
ent values of the nonlinear parameter a = 0, 0.2, 0.4 and 
0.6. For each plot, the solid line depicts the exact numer- 
ical results, whereas the dashed line corresponds to F^, 
Eq. ©J supplemented by a geometrical factor depending 
on 9. Away from the backward direction, the agreement 
between the exact numerical calculations and our theo- 
retical prediction for the background is clearly excellent. 
This is emphasized by the additional curve (long dashed 
line) plotted for a = 0.2 depicting the results obtained 
when neglecting the fluctuations of /(r) , for example re- 
placing (|ff(/)|2/) by |5((/))P (/) in Eq. ©. 

In the backward direction, constructive interference 
between reversed scattering paths results in the well- 
known coherent backscattering peak. As obvious from 
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Figure 2; Building blocks for the diagrammatic calculation of 
nonlinear coherent backscattering. Filled squares (with out- 
going solid arrows) denote the scattered field /, and open 
squares (with outgoing dashed arrows) the complex conju- 
gate /* . Incoming solid (dashed) arrows represent derivatives 
d/dE (d/dE*). Diagrams (a-d) contribute to the cooperon 
cross sections k and R, see Eq. ((Sjl, whereas diagrams (e,f) 
represent nonlinear propagation t* and r, see Eq. 0. 



Fig. [TJ the height of this peak is strongly affected by the 
nonhnearity. Nevertheless, we are perfectly able to in- 
corporate these interferences effects in our approach, see 
the horizontal dotted lines in Fig. [2 which depict the 
predicted total bistatic coefficient, F^ -I- Fr; see Eq. ([13]) 
below, in the exact backward direction These re- 

sults are obtained by a diagrammatic analysis, which we 
summarize in the following. 

As in the linear theory, we calculate the coherent 
backscattering effect by so-called 'crossed' or 'cooperon' 
diagrams describing pairs of reversed scattering 
paths. Hence, the individual scatterers are subject to 
two different incident probe fields E and E* , which rep- 
resent the two amplitudes propagating along the reversed 
paths. The response of a scatterer to these two weak 
probe fields is given by the derivative (P /{dEdE*). De- 
pending on whether the incident fields act on the dipole 
/ or its complex conjugate /*, we obtain the four terms 
represented in Fig. [2Ja-d). If we denote the sum of di- 
agram (a) -I- (c) by K, and (b) -|- (d) by k, the explicit 
expressions read: 
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If one of the incident fields originates from the coher- 
ent mode, d/dE is again replaced by l/E^, i.e. Kc = 
M{fdf*/dE*) /E, and =M{f*df/dE*) /E,. 

Concerning propagation between two scattering 
events, the refractive index, Eq. ([3]), remains unchanged 
for the reversed paths. In addition to that, however, we 
find two other contributions shown in Fig. [2)Je,f), which 
exist only as crossed diagrams. Here, diagram (e) repre- 



a) 

A 



b) 



A 



V 


t 






A 


A 



->n<r 
;a 

V 



Figure 3: Examples of forbidden (a) and allowed (b) combi- 
nations of the building blocks shown in Fig. [2] The forbidden 
combinations correspond to those forming a closed loop in the 
sense that two squares are connected by two arrows pointing 
in different directions. 



sents the expression: 
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and diagram (f) its complex conjugate (t*). For the co- 
herent mode, the above expression is again modified as 
follows: Tc = ~iN {d'^ f* l{dE*f) /{2kEc). The propa- 
gation of the thin line in Fig. [2l^e,f) is unaffected by the 
nonlinear event. 

The crossed transport equation is established by con- 
necting the building blocks shown in Fig. [2] with each 
other. As we have found, however, some combinations 
of diagrams represent unphysical processes which do not 
occur in the formal expansion of the solution of the non- 
linear wave equation as a multiple scattering series. An 
example is shown in Fig. [SJa). The problem with this 
diagram is that the fields radiated by /* and / mutu- 
ally depend on one another. Therefore, one cannot tell 
which one of the two events / or /* happens before the 
other one. This contradicts the multiple scattering series, 
where the individual scattering events occur one after the 
other. In order to avoid closed loops like the one shown 
in Fig. [3ja), we must ignore all combinations where one 
of the diagrams Fig.[2ljc,d) or (e) occurs after Fig.[2ljb,d) 
or (f) when following the solid arrow along the crossed 
path. 

We account for these forbidden diagrams by split- 
ting the transport equation into two parts, which we 
call Ci and C2. The first part, Ci, contains only dia- 
grams Fig. [DJa,c) and (e). As soon as one of the events 
Fig. [2fb,d,f) occurs, the crossed intensity changes from 
type Ci to type C2. The subsequent propagation of C2 
is then given by diagrams Fig. [51[a,b,f). Following these 
rules, we describe the propagation of C1.2 by transport 
equations similar to Eqs. (|4l5p : 



Cc{v) 



(10) 



Ci{r) = j^dv'P{Yy)(aCi+cj,c)j{Y') (11) 
C2(r) = j^dv'P{v,v')(a*C2+aCi+d,c}j(v')[12) 
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where the propagation kernel P(r, r') = 
exp(— |r — r'|/^)/(47r|r — r'|)^ is the same as in Eq. (O, 
and the cross sections cr result as follows: a — k + £Kt, 
a = k + £Kt* and, similarly, (Tc — + iKr,. and 
o'c — kc + (Kt* . Finally, the crossed bistatic coefficient 
reads: 

^ Jy ^e^'^("-"?)((a: + a:)Ci + (7:^2) (r) (13) 

For comparison with the background Fl, we define 
diffuson cross sections by writing K — a'^'^^Id + 
a^f' Ic^ such that Eq. ([5]) attains a form comparable to 
Eq. pT|) . Exploiting the Gaussian properties of the dif- 
fuse field, we find: a^"^^ = M {d{f f*) / dl) and a'f^ = 
Af {d{ff*)/dE) /E*. Thus, the decrease of the backscat- 
tering peak observed in Fig. [1] is traced back to the fact 
that the cooperon cross sections a and a decrease faster 
than the diffuson cross section a^'^\ Let us note that 
there also exist other models than Eq. ([7]) , where our the- 
ory predicts an increasing coherent backscattering cone. 
However, these models might suffer from speckle insta- 
bilities - a point which requires further investigations. 

To obtain the relatively simple form of Eqs. (|10m3p . 
we assume that the scattered intensity K{r) is approxi- 
mately constant on length scales comparable to the mean 
free path £, and we neglect some diagrams where the co- 
herent mode C'c is affected by a nonlinear event Tc. These 
approximations are expected to be well fulfilled in the 
case of large optical thickness b. In the numerical com- 
parison depicted in Fig. (TJ where b is not very large, we 
have used the exact version of Eqs. (|10m3p . which will 
be published elsewhere. 

As explained in the introduction, our theoretical 
scheme also applies to other types of nonlinear systems, 
like, for example, the case of linear scatterers embedded 
in a homogeneous nonlinear medium: 

AE{r) + (^e{r) + a\E{r)\^^E{r) = (14) 

with (5-correlated disorder e(r) corresponding to a (linear) 



mean free path Iq. Here, the dilute medium approxima- 
tion is vahd if k£o > 1 and {aI)'^k£o < 1. The latter 
condition is automatically fulfilled if we assume that we 
are in the stable regime, where Eq. (|14p has a unique so- 
lution. According to [2], this is the case (for a G R) if 
(a/)^6^(fc^o + b) < 1, with b the optical thickness. 

In this case, the diagrammatic method applies in the 
same way as described above. In particular, we obtain 
the following expressions for the cross sections: 

An 

a{r) = a,{r) = ^ [1 + tkioa (/c(r) + /^(r))] , (15) 
«o 

CT = CTc = -ATTika{Ic + Id), o-^'') = ai''^^ = Att/£o, and for 
the mean free paths n = (e) + a{Ic + Id) + i/{2k£o) and 
ric = {e) + a{Ic/2 + Id) + i/{2k£o). In the energy conserv- 
ing case a € R, it can be shown that C2 does not con- 
tribute to the real part of the backscattering coefficient 
Tc- Then, it follows from Eqs. ^ and that the 

nonlinearity introduces a phase difference A0 — Mki^al 
between reversed paths undergoing M linear scattering 
events. Since (M) oc 6, we predict a significant reduction 
of the coherent backscattering peak if bkl^al ~ 1 (which 
is still inside the stable regime if kto is large). 

In summary, we have extended the usual diagrammatic 
approach to take into account nonlinear effects for the co- 
herent transport in disordered systems beyond the per- 
turbative regime. The excellent agreement with direct 
numerical simulations emphasizes the validity of our ap- 
proach. It readily applies for many different nonlinear 
wave equations. Eq. (fT4|) . for example, is mathemati- 
cally equivalent to the Gross-Pitaeskii equation describ- 
ing nonlinear propagation of matter waves in random po- 
tentials. In the latter case, this method will allow us 
to describe the localization properties of the mean field. 
Extending the present approach within the Bogolioubov 
framework, it will be possible to understand how these 
localization properties are affected by the non-condensed 
fraction of the atoms. 
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